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Abstract 

We  obtain  analytical  solutions  for  the  perturbed  shock 
paths  induced  by  time-varying  random  motions  of 
a  piston  moving  inside  an  adiabatic  tube  of  con¬ 
stant  area.  The  variance  of  the  shock  location  grows 
quadratically  with  time  for  early  times  and  switches  to 
linear  growth  for  longer  times.  The  analytical  results 
are  confirmed  by  stochastic  numerical  simulations,  and 
deviations  for  large  random  piston  motions  are  estab¬ 
lished. 

The  English  physicist  James  Joule  was  perhaps  the 
first  to  use  the  concept  of  a  moving  piston  in  order 
to  demonstrate  the  mechanical  equivalent  of  heat  in 
his  pioneering  studies,  almost  two  centuries  ago.  The 
moving  piston  has  also  been  used  extensively  in  funda¬ 
mental  studies  of  fluid  mechanics  and  shock  disconti¬ 
nuities  in  the  last  century,  and  this  now  classical  prob¬ 
lem  has  been  solved  analytically  in  one-  and  also  higher 
space  dimensions  (1,  2).  It  is  well  known  that  a  shock 
wave  propagating  into  a  stationary  fluid  sets  it  into 
motion  and  raises  its  pressure,  temperature  and  den¬ 
sity.  This  situation  can  be  physically  realized  by  a 
planar,  cylindrical  or  spherical  piston  moving  at  spec- 
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ified  speed  into  a  stagnant  fluid.  In  gasdynamics,  in 
particular,  in  the  context  of  normal  shock  waves,  the 
one-dimensional  classical  problem  describes  a  piston 
moving  at  constant  speed  in  a  tube  of  constant  area 
and  adiabatic  walls;  the  shock  wave  is  created  ahead 
of  the  piston.  Solutions  of  this  flow  problem  with  arbi¬ 
trary  piston  speeds  are  difficult  to  obtain  even  for  the 
one-dimensional  case,  although  recently  approximate 
analytical  solutions  have  been  obtained  for  accelerat¬ 
ing  and  decelerating  pistons  valid  only  for  short  times, 
e.g.  see  (3). 

In  the  current  work  we  revisit  the  one-dimensional 
piston  problem  within  the  stochastic  framework,  i.e., 
we  allow  for  random  piston  motions  which  may  be 
changing  in  time.  In  particular,  we  superimpose  small 
random  velocity  fluctuations  to  the  piston  velocity  and 
aim  to  obtain  analytical  solutions  of  the  stochastic  flow 
response.  Within  the  context  of  small  random  fluctu¬ 
ations,  we  assume  that  the  same  thermodynamic  con¬ 
ditions  are  valid  as  in  the  classical  problem,  i.e.,  that 
an  isentropic  region  exists  between  the  piston  surface 
and  the  shock  wave.  This  assumption  is  justified  by 
the  theory  of  weak  shocks  although  at  the  microscopic 
level  more  complex  processes  may  take  place.  For  ex¬ 
ample,  it  has  been  reported  in  (4)  that  a  wall,  which 
is  adiabatic  when  rigidly  fixed,  may  become  conduct¬ 
ing  when  it  is  allowed  to  have  a  stochastic  motion  in¬ 
dependently  of  the  value  of  its  macroscopic  velocity. 
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However,  in  the  macroscopic  models  we  develop  here 
we  will  assume  that  such  effects  are  negligible  and  thus 
all  surfaces  remain  adiabatic. 

In  the  first  part  of  the  paper,  we  employ  stochas¬ 
tic  perturbation  analysis  to  obtain  closed-form  ana¬ 
lytical  formulas  for  the  perturbed  shock  paths.  The 
random  piston  motion  is  modeled  as  a  stochastic  pro¬ 
cess  following  a  Markov  chain  corresponding  to  various 
values  of  correlation  length.  The  main  physical  find¬ 
ing  extracted  from  the  analytical  solution  is  that  the 
variance  of  the  location  of  the  perturbed  shock  grows 
quadratically  with  time  at  early  times  but  it  switches 
to  linear  growth  at  later  times.  In  the  second  part 
of  the  paper,  we  perform  high-resolution  stochastic 
simulations  using  a  standard  Monte  Carlo  approach 
and  also  using  the  polynomial  chaos  method  based  on 
Wiener-Hermite  expansions.  The  objective  is  to  con¬ 
firm  the  results  of  perturbation  analysis  and  determine 
their  validity  range  using  numerical  solutions  of  the 
full  nonlinear  Euler  equations  subject  to  stochastic  in¬ 
puts.  More  generally,  the  stochastic  piston  problem  we 
have  defined  here  serves  as  a  strict  testbed  for  rigorous 
evaluation  of  numerical  stochastic  solvers,  and  to  this 
end,  we  have  compared  the  performance  of  polynomial 
chaos  against  the  Monte  Carlo  approach.  The  results 
depend  critically  on  the  specific  value  of  correlation 
length  as  well  as  on  the  length  of  time  integration.  At 
early  times  and/or  large  values  of  correlation  length 
the  polynomial  chaos  method  outperforms  (often  by 
orders  of  magnitude)  the  Monte  Carlo  approach,  how¬ 
ever  it  is  not  as  effective  in  other  cases. 


determined  in  terms  of  the  piston  speed  through  the 
conservation  of  mass,  momentum  and  energy  (5).  For 
perfect  gas  with  constant  specific  heats,  these  relations 
are: 
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where  7  =  cp/cv  is  the  ratio  of  specific  heats,  and 
Co  is  the  local  sound  speed  ahead  of  the  shock.  The 
sound  speed  behind  the  shock  can  be  obtained  from 
Eqs.  la  and  lb: 

C2  =  (C20+7SUp)(1-Up/S).  (2) 


Figure  1:  Sketch  of  piston-driven  shock  tube  with  ran¬ 
dom  piston  motion. 


1  Stochastic  Perturbation  Anal¬ 
ysis 

We  consider  a  piston  having  a  constant  velocity,  Up, 
moving  into  a  straight  tube  filled  with  a  homogeneous 
gas  at  rest.  A  shock  wave  will  be  generated  ahead  of 
the  piston.  A  sketch  of  the  piston-driven  shock  tube 
with  random  piston  motion  superimposed  is  shown  in 
Fig.  1.  Given  the  state  ahead  of  the  shock,  the  speed 
of  the  shock,  S,  and  the  thermodynamic  states  of  the 
gas  behind  the  shock  (i.e.,  ahead  of  the  piston)  are 


In  the  following  we  will  normalize  all  velocities  with 
Co,  and  thus  Co  =  1.  We  now  define  the  stochastic  mo¬ 
tion  of  the  piston  by  superimposing  a  small  stochastic 
component  to  the  steady  speed  of  the  piston,  i.e., 

up(t)  =  Up[l  +  eV(t,u>)],  (3) 

where  the  amplitude  e  is  small,  0  <  e  <C  1.  Our  objec¬ 
tive  is  to  find  how  do  the  perturbed  shock  paths  due 
to  the  random  piston  motion  deviate  from  the  unper¬ 
turbed  ones;  the  latter  are  given  by 

X(t)  =  S  ■  t.  (4) 
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Under  the  small  amplitude  assumption,  the  flow  field 
induced  by  this  perturbation  can  be  obtained  based 
on  the  assumption  that  the  propagation  speed  in  the 
region  behind  the  shock  and  ahead  of  the  piston  can  be 
identified  as  the  propagation  speed  of  the  unperturbed 
flow  quantities,  i.e.,  UP±C,  where  C  is  the  unperturbed 
sound  speed  behind  the  steadily  moving  shock  as  given 
by  Eq.  2. 


Figure  2:  A  sketch  of  shock  paths  induced  by  random 
piston  motions. 

To  proceed  we  consider  the  perturbed  Riemann’s 
invariants  and  evaluate  them  at  the  shock  from 

•  _l_  ^  .tl  ,  „ S  +  S'Up 

(l±k)vs  =  j±  =  vp± - -ap,  with  k  =  C  - — - — 

7  ~  1  1  +  jSUp 

(5) 

where  S'  =  and  vp,  ap  and  j±  are,  correspond¬ 
ingly,  the  perturbed  piston  velocity,  the  perturbed  lo¬ 
cal  sound  speed,  and  the  perturbed  Riemann’s  invari¬ 
ants.  These  invariants  are  constant  along  the  unper¬ 
turbed  (straight)  characteristic  lines.  For  more  details 
on  this  derivation  and  on  the  assumptions,  the  inter¬ 
ested  reader  is  referred  to  (6).  Fig.  2  shows  a  sketch 
of  the  shock  paths  induced  by  random  piston  motions. 
Specifically,  the  distorted  lines  show  an  instantaneous 
realization  of  the  piston  path  and  shock  path.  They 
are  distorted  due  to  induced  reflections  as  sketched  in 
the  plot  via  the  characteristic  lines.  In  the  sketch, 


the  steady  and  perturbed  piston  paths  are  denoted 
by  Upt  and  Upt  +  e  r/(t)  while  those  for  the  shock 
paths  by  St  and  St  +  e  £(t).  Also,  Vpfon+i)  and 
vs(t2n)  are  on  the  forward  characteristic  %  =  Up  +  C; 
vs{t2n+2)  while  vp(t2n+i)  is  on  the  backward  charac¬ 
teristic  =  Up  —  C.  Thus,  we  have  through  the  use 
of  the  Riemann  invariants  in  Eq.  (5),  that: 

2 

(1  +  k)vs(t2n)  =  Vp(t2n+i )  H - Tap(t2n+i)  (6a) 

7-1 

2 

(1  -  k)vs(t2n+ 2)  =  vp(t2n+i) - -ap(t2„+i).  (6b) 

7-1 

Adding  Eqs.  6a  and  6b  to  eliminate  o,p(t'2n+i),  we 
obtain  the  following  recurrence  formula: 

vs(t2n)  =  qvp(t2 „+i)  -  rvs(t2n+ 2),  n  =  0, . . . ,  N, . . . 

(7) 

where 

2  ,  1  -  k 

q  =  - - -  and  r  =  - - - . 

1  +  k  1  +  k 

Eq.  7  defines  a  recursive  relationship  between  the  ve¬ 
locities  at  the  shock  and  the  perturbation  of  the  piston 
motion  vp{t).  Starting  at  time  to  and  iterating  up  to 

N,  we  obtain  from  Eq.  7  a  set  of  (fV  +  1)  terms.  Elimi¬ 
nating  vs(t2),  vs(t4) .  ■  ■  vs(t2N)  from  this  set  we  obtain 

N 

vs(t)  =  q^2,(-r)nvp{t2n+i)  +  {-r)N+1vs(t2N+ 2)- 

n— 0 

(8) 

If  the  perturbation  of  the  piston  starts  at  time  ts  > 

O,  the  zigzag  path  of  the  characteristics  coming  down 
to  the  origin  will  end  on  the  piston  path;  therefore,  vs 
in  the  last  term  in  Eq.  8  is  zero.  On  the  other  hand, 
if  ts  =  0,  the  zigzag  path  will  zigzag  indefinitely  to 
approach  t  =  0  i.e.,  N  — >  00.  Since  r  is  always  less 
then  unity,  the  last  term  of  Eq.  8  will  approach  zero 
for  any  finite  value  of  vs(t0 0).  Therefore,  one  can  drop 
the  last  term  in  Eq.  8  to  obtain: 

N 

vs(t)  =  qj^(-r)nvp(t2n+ 1),  (9) 

n— 0 

where  N  =  00  if  ts  =  0,  i.e.,  the  perturbation  of  the 
piston  starts  at  t  =  0  or  it  is  determined  by  the  last 
non-zero  value  of  vp(t2N+i)- 
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To  make  effective  use  of  this  relation  we  need  to 
obtain  the  relationship  of  the  shock  locations  at  t2n, 
t2n+i  and  t2n+2-  To  this  end,  let  us  denote  the 
shock  locations  at  time  t2n  and  f2n+2  by  Xs(t2n)  and 
Xs(t2n+2),  correspondingly,  and  the  piston  location  at 
t2n+i  by  Xp(t2n+i)-  Assuming  the  characteristics  are 
approximated  by  straight  lines  with  their  slopes  given 
by  Up  ±  C,  we  have 

Xs(t2n')  Ap(f2n+l)  —  ( Up  -T  C)  (t2n  t2n+l) 

Xp(t2n+l)  ~  Xs(t2n+2)  =  (Up  —  C) 

(t-2n+l  —  t2n+2)-  (10) 

Also,  defining  the  perturbed  path  of  the  shock  and  the 
piston  by  £(t)  and  r?(t),  we  can  express 


Eq.  8  then  reads  as 

N 

vs(t)  =  q^2(-r)nvp(aPnt),  (12) 

n= 0 

and  the  shock  speed  is  then  obtained  as 

7  Q 

S(UP  +  vs(t))  =  S(Up)  +  -jjj-vs(t)  =  S  +  S'vs(t), 
with  the  shock  path  governed  by 


Using  Eq.  11  we  have 


Xs(t2n)  =  S  ■  t,2n  +  U;(t2n) , 
Xs(t2n+2)  =  S  '  t.2n+2  +  eti(t2n+2) , 
Xp(t2n+l)  =  Upt2n+1  +  CT/fen+l)- 


Finally,  substituting  these  into  Eq.  10,  we  obtain 

t2n+l  =  cd2n  +  [vfan+l)  ~  ^(^2n)] 

^2n+2  =  Pt2n  +  tlfivifan+l)  ~  ^(^2n) 

-  C(^2n+2)],  (11) 


where 


a 

P 
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C  +  Up -s 
c 

C  +  Up  -  s 
C  +  S  -Up 

1 

c  +  s-up 


The  inequalities  above  are  due  to  C  >  1  and  S  >  Up. 

If  we  drop  the  term  containing  e  in  Eq.  11,  which 
is  consistent  with  the  small  disturbance  assumption, 
the  recurrent  relationship  is  very  much  simplified  and 
closed  form  solutions  can  be  obtained.  (We  will  retain 
this  term  in  the  next  subsection  below).  With  this 
simplification  we  have  that: 


t2n+l  —  crfgn  —  Ot(3t2n-2  —  Ot(32t2n-l  ■  ■  ■  —  a/3Uto 

=  a/3nt. 


ef =QS'Y/(-r)nVp(apnt ),  (13) 

71=0 

and  taking  £(0)  =  0,  we  obtain 

O/  00  nt 

m  =  —  £(-r)"  /  dhvpia^h).  (14) 

6  n—0  1'0 

As  a  check  of  Eq.  12,  we  first  consider  the  simple 
problem  of  a  piston  with  its  velocity  subject  to  a  small 
constant  perturbation  vp  starting  at  t  =  0  (e.g.,  a  step 
function  of  size  vp).  According  to  Eq.  13,  we  have 

OO  OO 

vs(t)  =  g  (— r)nUp(£2rt+i)  =  qvp  ^2  ~r"  =  VP 

71=0  71=0 

i.e.,  the  velocity  behind  the  shock  will  be  Up  +  vp  for 
all  t  >  0.  The  shock  speed  S  will  change  to  S(Up  +  vp), 
instead  of  S(UP). 

We  now  consider  vp  to  be  a  random  process  with 
zero  mean  and  the  following  covariance 

i ’p(t)  =  eUpV(t,io) 

{V(t,w))=0 

(V(ti,u),V(t2,uj))  =  cr '  A  2' 


(15a) 

(15b) 
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where  A  is  the  correlation  time.  The  above  covariance 
kernel  describes  a  Markov  random  process  in  time. 
The  larger  the  value  of  the  correlation  time  A  is  the 
closer  the  random  motion  approaches  a  fully-correlated 
process  -  we  refer  to  this  as  a  random  variable  case.  On 
the  other  hand,  the  smaller  the  value  of  the  correlation 
time  A  is  the  closer  the  motion  resembles  white  noise. 

Substituting  Eq.  15a  into  14,  we  obtain 


N  ft 

m=qS'Upy^{~r)n  dtiV(a(3nti,uj). 

n—0 


Because  of  Eq.  15b,  we  have 


m)  =  o 

oo  oo 

<m  =  M2EE(-r)m+" 

71=0  771=0 

f  dt i  f  df2e-3l/3n‘i-/3m^  (16) 

Jo  Jo 

The  double  summation  in  Eq.  16  can  be  split  into 
three  parts:  1.  sum  of  all  diagonal  terms,  2.  sum  of 
all  the  terms  above  the  diagonal,  and  3.  the  terms 
described  in  detail  below.  It  is  easy  to  see,  that  the 
last  two  sums  are  equal.  Thus,  we  have 
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<r(t)>  =  (UPqS') 


OO  71—1 


71+771 


dt\ 


71=1  771  =  0 


[  dt2e~^ntl~0mt21 

Jo 

oo  »t  ot 

+  V(r2n)  /  dti  /  dt2e 
Jo  Jo 


-s+|q-t2i 


Both  integrals  in  the  above  equation  can  be  integrated 
explicitly  to  give 


(£2(t)) 


{UpqS'A/af 


OO  71—1 

2EE(-r)"+m^M+ 

71=1  771=0 


EE  7'2ninAT) 

71=0 


(17) 


where  r  =  at  I  A,  and 

e~0r"T  +  e~^T  -  1  -  e-(/3m-/3n)r 

where  m  <  n.  For  r  <C  1  it  is  easy  to  show  that 
In,m  =  t2 .  The  summations  in  Eq.  17  can  be  per¬ 
formed  explicitly  to  obtain 

<?2(r))  «  (UpqS'A/af  for  T  «  i.  (i8) 


.w  = 


1 


ZT 


am  an+m 


For  r  1,  we  can  neglect  the  exponential  terms  in 
equation  (1),  and  thus 


In,m(T )  — 


2  r 

/3m 

2  r 

+ 


1 

^71+771 


for  m  <  n 


In 


for  m  = 


These  expressions  for  /„jm  can  again  be  summed  ana¬ 
lytically  on  the  right-hand-side  of  Eq.  17,  to  obtain 


<£2(r)> 


( UpqS'A/a )2 


2r(l  —  r) 

(1  -  r2//3)(  1  +  r) 


2 

(1  —  r2//32)(l  +  r//3) 


for  r  1 


(19) 


For  arbitrary  values  of  r,  we  calculate  the  quan¬ 
tities  in  the  square  brackets  in  Eq.  17  numerically. 
Because  of  the  smallness  of  the  values  of  r  and  /3,  the 
series  converges  fast.  We  plot  in  Fig.  3  the  quantity 
(?2(r)) /(UpqS' A/ a)2  as  a  function  of  r  given  by  Eq. 
17  with  Up  =  1.25  i.e.,  corresponding  to  Mach  number 
of  the  shock  M  —  2.  The  asymptotic  formula  for  small 
and  large  r  given  in  Eqs.  18  and  19  are  also  included 
in  the  plot.  We  observe  a  qualitative  change  in  the 
stochastic  response  versus  time.  At  early  times,  the 
location  of  the  path  scales  linearly  with  time  whereas 
at  later  times  it  scales  with  square  root  of  time  (note 
that  the  variance  <  £  >2~  (Length)2).  This  interest¬ 
ing  result  is  consistent  with  physical  intuition  suggest¬ 
ing  that  at  early  times  convective  motions  dominate 
while  at  longer  times  the  diffusion  process  takes  over. 
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Figure  3:  Normalized  variance  of  perturbed  shock 
paths.  Solid  line:  perturbation  analysis  results,  see 
equation  (17);  Dashed  line:  early-time  asymptotic  re¬ 
sults  from  Eq.  18;  Dash-Dotted  line:  late-time  asymp¬ 
totic  results  from  Eq.  19. 


n,  we  have  the  following  sequence: 

t'2n  =  Pt2(n-1)  +  ej[2r](t2n-l)  -  €(t2(„- 1)) 

-£(*2n)] 

^2(n—l)  =  /^2(n-2)  +  ^[^V(t2n-3)  ~  £(^2(n-2)) 

-£(*2(n-l))] 

t2(n-2)  =  Pt2(n-3)  +  e7[2??(i2n-5)  —  €(t2(n-3)) 

^(^2(71.— 2)  )] 

t2  =  /3t  +  e'/[2ri(t1)  -  £(t)  -  €(t2)\-  (20) 

Solving  for  t2n  from  the  above  equations,  we  obtain: 

n 

t2n  =(3nt  +  eq[E  -  (1  +  /3)  (21) 

i= i 

n—1 

1=1 

Thus,  t2n+i  can  be  expressed  as: 

t2n+ i  =  (3n[at  -  +  ^[2/?E^(^-i)/3n-J  + 

i=i 


V(t2n+l)\  ~  E  (22) 

Within  the  context  of  small  amplitude  random  mo-  J_1 

tions  of  the  piston,  we  have  neglected  the  last  term  in  ^he  sh°ck  path  is  governed  by 
Eq.  11  involving  e.  This  allowed  us  to  obtain  closed-  (]^ 

form  analytical  solutions,  as  we  explained  in  the  pre-  =  iS  vp(J 2n+i), 

vious  section.  Now  we  revisit  this  approximation  and  n~° 

retain  that  term,  so  we  employ  the  following  recurrence  and  taking  £(0)  =  0,  we  have 
formulas: 

n/  00  ft 

m  =  —  E(-r)"/  ^i^(wi)-  (23) 
6  n= 0  "'° 

Considering  Eqs.  15b  and  23  together,  we  obtain 
^2n+l  =  Oit2n  +  fj[r)(t2n+l)  ~  £(t2n)}  (£(t))  =  0 

t'2n+2  =  Pt2n  +  £7[277(<2n+l)  —  £(t2n)  ~  £(t2n+2)\ 

OO  OO  pt 

<m>  =  M2EEHm+"  d*i 

n=0  m—0  •'° 

where  a,  (3  and  7  are  given  in  Eq.  (12).  For  a  general 


/  dt2e_  *  l*1'2'i+1_*2’2m+1I 

Jo 


6 


To  compute  the  variance  of  the  induced  shock  path 
(£2(f))  we  need  to  compute  t^n+i-  However,  Eq.  22 
shows  that  to  compute  t2n+i  we  have  to  know  the 
shock  path  £(t)  at  all  the  previous  reflection  times 
t.2j+i  and  t2j-  To  this  end,  we  solve  Eqs.  21  22,  23 
and  24  numerically,  by  employing  an  iteration  method 
and  setting  t.2j+ i  =  afPt  and  t2j  =  /3H  as  an  initial 
approximation . 

In  the  following  section,  we  perform  stochastic  nu¬ 
merical  simulations  to  confirm  our  findings  and  estab¬ 
lish  limitations  of  the  stochastic  perturbation  analysis 
presented  in  this  section. 


2  Stochastic  Simulations 

We  perform  two  types  of  stochastic  simulations  to 
verify  the  results  of  the  previous  section,  following 
a  Monte  Carlo  approach  and  a  polynomial  chaos  ap¬ 
proach.  We  employ  the  full  nonlinear  Euler  equations 
with  the  extra  complication  that  there  is  an  unsteady 
stochastic  boundary,  namely  the  piston  position.  To 
this  end,  a  boundary-fitted  coordinate  approach  is  em¬ 
ployed  to  transform  the  equations  into  a  stationary 
domain.  The  transformed  Euler  equations  contain 
stochastic  source  terms  proportional  to  —pdup/dt  and 
—pvdup/dt  in  the  momentum  and  energy  equations, 
respectively. 

In  Monte  Carlo  simulations,  we  use  a  Markov  chain 
in  time  to  represent  the  stochastic  input.  In  polyno¬ 
mial  chaos  simulations,  the  representation  of  stochastic 
inputs  is  expressed  by  a  Karhunen-Loeve  decomposi¬ 
tion;  see  references  (7,  8).  Specifically,  we  consider 
different  representations  of  the  stochastic  inputs  cor¬ 
responding  to  the  covariance  kernel 

<  V(ii,w)V(i2,u;)  >=  (25) 

where  A  is  the  correlation  length.  A  corresponding 
Markov  chain  is  employed  to  represent  discretely  the 


exponential  kernel  as  follows: 

Vo  =  Co 
Vi=  C%  +  /C  i 

Vi+\  =  CVi  +  ,/Ci+i 

where 

C=  e^1  and  /=  \]\  -  C 2. 

In  the  Monte  Carlo  simulation,  a  random  piston  ve¬ 
locity  Up  =  Up(l  +  eVi(t,u>))  is  selected  from  the  above 
Markov  chain  as  a  stochastic  input  at  each  time  step 
t,.  In  the  polynomial  chaos  representation  we  em¬ 
ploy  Wiener-Hermite  expansions  for  all  conservative 
and  derived  stochastic  variables  of  the  form 

M 

X{w)=YiZj*Mu))>  (26) 

3=0 

where  the  basis  {4>j}  is  formed  from  the  Hermite 
orthogonal  polynomials  of  degree  p.  Here  £(w)  is 
a  Gaussian  variable  of  dimension  N  and  M  is  the 
total  number  of  deterministic  coefficients  Xj,  where 
M  +  1  =  ( N  +  p)\/(N\p\).  We  employ  the  fifth-order 
WENO  method  in  space  in  order  to  capture  the  shock 
location  accurately  and  the  third-order  TVD  R.unge- 
Kutta  method  in  time;  see  details  in  (9). 

We  now  present  some  results  for  the  following  con¬ 
ditions:  Behind  the  shock  we  impose  a  steady  piston 
velocity:  Up  =  1.25  (normalized  by  the  sound  speed 
ahead  of  the  shock),  i.e.,  corresponding  to  Mach  num¬ 
ber  of  the  shock  M  =  2.  Ahead  of  the  shock  the  sound 
speed  is  C0  =  1  and  the  pressure  is  P  =  1.  We  inves¬ 
tigate  the  stochastic  response  for  various  values  of  the 
correlation  length  A  and  of  the  amplitute  of  the  ran¬ 
dom  piston  motion  e. 

In  Fig.  4,  we  plot  the  variance  of  the  perturbed 
shock  paths  induced  by  small  random  piston  motions 
corresponding  to  amplitude  e  =  0.01  and  correlation 
lengths  A  =  0.5, 1, 2  and  10  obtained  from  Monte 
Carlo  simulations  (2,000  runs).  There  is  good  agree¬ 
ment  of  the  Monte  Carlo  solutions  with  the  analytical 
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Figure  4:  Solid  line:  Variance  of  the  perturbed  shock 
paths  for  e  =  0.01  and  A  =  0.5, 1,2  and  10  obtained 
from  Monte  Carlo  simulations;  Dashed  line:  results 
from  perturbation  analysis;  Dash-Dotted  line:  results 
from  equation  (19)  for  later  time;  Dash-Dot-Dotted 
line:  results  from  equation  (18)  for  early  time. 

solutions.  In  Fig.  5,  we  plot  the  variance  of  the  per¬ 
turbed  shock  paths  induced  by  random  piston  motions 
corresponding  to  correlation  length  A  —  1  and  am¬ 
plitudes  e  =  0.01,0.1,0.2,0.3  and  0.5  obtained  from 
Monte  Carlo  simulations  (3,000  runs).  For  small  am¬ 
plitudes  e  =  0.01,0.1  and  0.2,  good  agreement  is  ob¬ 
served  between  Monte  Carlo  simulations  and  analyt¬ 
ical  solution.  However,  for  larger  amplitude,  such  as 
e  =  0.3  and  0.5,  the  stochastic  simulation  deviates 
from  the  analytical  solution.  We  will  examine  this 
discrepancy  in  some  more  detail  below  but  first  we 
present  results  from  the  polynomial  chaos  simulations. 

Fig.  6  shows  results  from  polynomial  chaos  simula¬ 
tions  corresponding  to  piston  motions  described  by  a 
random  variable,  i.e.  a  fully-correlated  stochastic  pro¬ 
cess  whereby  A  — >  oo.  The  polynomial  chaos  simula¬ 
tions  match  quite  closely  the  exact  analytical  solutions 
even  over  a  more  than  two-orders  of  magnitude  change 
in  the  value  of  the  variance.  This  verifies  the  conver¬ 


Figure  5:  Dotted  line:  Variance  of  the  perturbed  shock 
paths  for  A  =  1  and  e  =  0.01,0.1,0.2,0.3  and  0.5  ob¬ 
tained  from  Monte  Carlo  simulations;  Solid  line:  re¬ 
sults  from  perturbation  analysis;  Dashed  line:  results 
from  equation  (19)  for  later  time;  Dash-Dotted  line: 
results  from  equation  (18)  for  early  time. 


gence  of  Hermite-chaos  for  this  case.  Fig.  7  shows 
results  from  polynomial  chaos  simulations  correspond¬ 
ing  to  piston  motions  described  by  a  random  process 
with  amplitude  e  =  0.01  and  correlation  time  A  =  1. 
In  the  polynomial  chaos  simulations,  the  number  of 
stochastic  dimensions  of  random  input  is  changed  from 
N  =  3,6, 50  to  100  (N  is  also  the  number  of  Karhunen- 
Loeve  modes  for  representing  the  stochastic  piston  mo¬ 
tion).  By  increasing  the  dimensions  of  random  input, 
the  polynomial  chaos  simulations  agree  better  with  the 
analytical  solution  longer.  However,  there  is  a  finite  er¬ 
ror  after  long  time  integration  unlike  the  Monte  Carlo 
simulations. 

We  now  re-examine  what  is  the  effect  of  neglecting 
the  second  term  in  Eq.  (11),  which  we  included  in  the 
refined  perturbation  analysis  of  the  previous  section. 
In  Fig.  8,  we  compare  the  variance  of  the  perturbed 
shock  paths  with  large  Up  +  vp{t)  random  piston  mo- 


Figure  6:  Variance  of  the  perturbed  shock  paths  for  a 
random  variable  (fully-correlated  kernel,  A  — >  oo)  with 
amplitude  e  =  0.01.  Dash-Dotted  line:  analytical  solu¬ 
tion  from  perturbation  analysis;  Solid  line:  numerical 
results  from  polynomial  chaos  simulations. 


Figure  7:  Variance  of  the  perturbed  shock  paths  for 
e  =  0.01  and  A  =  1  obtained  from  polynomial  chaos 
simulations  with  stochastic  dimensions  N  =  3,6,50 
and  100. 


tions  obtained  from  Monte  Carlo  simulations,  analyt¬ 
ical  solutions  from  perturbation  analysis,  and  analyti¬ 
cal  results  obtained  including  the  corrections  for  larger 
random  piston  motions.  Indeed,  significant  improve¬ 
ment  in  the  semi-analytical  results  is  evident  compared 
to  Monte  Carlo  simulations. 


Figure  8:  Thick-Solid  line  and  Thick-Dashed  line: 
Variance  of  the  perturbed  shock  shock  paths  corre¬ 
sponding  to  correlation  length  A  =  1,  amplitudes 
e  =  0.1  and  0.3  obtained  from  larger  random  piston 
motion  modeling;  Dotted  line:  results  from  Monte 
Carlo  simulations;  Solid  line:  results  from  perturba¬ 
tion  analysis;  Dashed  line:  results  from  equation  (19) 
for  later  time;  Dash-Dotted  line:  results  from  equation 
(18)  for  early  time. 


3  Summary 

The  stochastic  piston  problem  is  a  re-formulation, 
within  the  stochastic  framework,  of  a  classical  aerody¬ 
namics  problem  that  studies  how  small  random  piston 
motions  affect  shock  paths.  We  have  developed  an  an¬ 
alytical  solution  for  the  linearized  Euler  equations  for 
the  stochastic  piston  problem.  Specifically,  Eqs.  (17), 
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(18)  and  (19)  represent  the  main  analytical  results  of 
this  paper.  The  first  equation  gives  the  full  analytical 
expression  whereas  the  last  two  give  asymptotic  results 
for  early  times  and  longer  times,  respectively.  They  re¬ 
veal  that  the  variance  of  the  location  of  the  perturbed 
shock  paths  grows  initially  quadratically  with  time  and 
switches  to  linear  dependence  for  longer  times.  The 
stochastic  numerical  simulations  presented  in  the  pa¬ 
per  confirm  the  results  and  show  good  agreement  with 
the  analytical  solution  for  up  to  20%  amplitudes  of  the 
random  piston  motion  compared  to  the  mean  steady 
motion. 
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